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Background: The current understanding of the earthquake interevent times distribution (ITD) in 
terms of seismological laws or fundamental physical principles is incomplete. The Weibull distri- 
bution is used to model the earthquake ITD. 

Purpose: To better understand the earthquake ITD in terms of fracture mechanics. 

Method: We link the earthquake ITD on single faults with the Earth's crustal shear strength dis- 
tribution, for which we use the Weibull model, by means of a phenomenological, stochastic stick - 
slip model. We generalize the model for fault systems. 

Results: We obtain Weibull ITD for power-law stress accumulation, i.e., a(t) = a t 13 , where /3 > 
for single faults or fault systems with homogeneous strength statistics. We also show that log- 
arithmic stress accumulation leads to the log-Weibull ITD. For the Weibull ITD, we prove that (i) 
m = P m s , where m and m s are, respectively, the interevent times and crustal shear strength Weibull 
moduli and (ii) the time scale t s = (Sg/a) 1 ^ where S s is the scale of crustal shear strength. We in- 
vestigate deviations of the ITD tails from the Weibull due to sampling bias, magnitude selection, and 
non-homogeneous strength parameters. Assuming the Gutenberg - Richter law and independence of 
m on the magnitude threshold, Ml )C; we deduce that t s oc e~ pMML ' c , where pu £ [1-15, 3.45] for 
seismically active regions. We demonstrate that a microearthquake sequence conforms reasonably 
well to the Weibull model. 

Conclusions: The stochastic stick - slip model justifies the Weibull ITD for single faults and for 
homogeneous fault systems, while it suggests mixtures of Weibull distributions for heterogeneous 
fault systems. Non-universal deviations from Weibull statistics are possible, even for single faults, 
due to magnitude thresholds and non-uniform parameter values. 
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I. INTRODUCTION 

Tectonic earthquakes are the result of the stick - slip motion of plates within the Earth's crust. 
This motion can be viewed in the framework of driven dissipative systems. Earthquakes are com- 
plex processes that involve nonlinearities, stochasticity, and multiple spatiotemporal scales. Earth- 
quakes originate on individual faults, which can be viewed as the fundamental structural units for 
earthquake processes. Faults are geological structures that span a variety of scales, from a few cen- 
timeters up to thousands of kilometers for plate boundary faults, and accommodate earthquakes 
that span up to eight orders of magnitude in size [HI El. Neighboring faults interact with each other 
leading to space-time organization of earthquakes within fault systems, e.g. Q. The statistical 
laws of seismicity, e.g., the Gutenberg - Richter scaling law, are valid for systems of faults but not 
necessarily for single faults flU. 

One of the most interesting problems in statistical seismology and statistical physics is the 
probability distribution of the times between successive earthquake events that exceed a specified 
magnitude threshold, i.e., the so-called interevent times distribution. The terms return intervals, 
waiting times and recurrence intervals are also used instead of interevent times. A subtle but im- 
portant distinction can be made between recurrence intervals, which refer to consecutive events 
that take place on the same fault, and interocurrence intervals, which focus on all faults in a spec- 
ified region [5|. The statistical properties of recurrence intervals are more difficult to estimate, be- 
cause less information is available for individual faults. The distinction is, nevertheless, conceptu- 
ally important, since recurrence intervals characterize the one-body, i.e., the single-fault, problem, 
while interoccurrence intervals are associated with the activity of the many-body system 0. In the 
following we use the term interevent times in both cases, but we distinguish between models that 
apply to the one-body versus the many-body problems. 

Knowledge of the interevent times distribution is important for the assessment of earthquake 
hazard in seismically active regions. From a practical perspective, inferring interevent times statis- 
tics for very large earthquakes based on the statistics of smaller-magnitude and higher-frequency 
events is highly desirable. In seismology, earthquakes are categorized as foreshocks, main shocks 
and aftershocks. Aftershocks are assumed to be triggered by dynamic stress changes induced by 
the main event, such as stress redistribution, fluid flow triggering, etc., that are not directly linked 
to the tectonic motion of the plates. The classification of earthquakes as main events, foreshocks 
and aftershocks (declustering), is based on heuristic arguments. The result of declustering anal- 
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ysis strongly depends on the method used @. Given this ambiguity and the fact that all seismic 
events in a fault system can be viewed as the result of the same dynamic process, we do not dis- 
tinguish between different types of events; this approach is commonly used in statistical physics 
studies GHJH. 

A. Notation 

We briefly comment on some notation and abbreviations used in the following: the probability 
density function of a random variable X(t) is denoted by fx(%) and the cumulative probability 
function by Fx{x) = P(X < x). The strength of the Earth's crust is denoted by the random vari- 
able S, which is assumed to follow Weibull statistics with shape parameter m s and scale parameter 
S s . The local magnitude (Richter scale) of an earthquake event will be denoted by M L , and the 
magnitude threshold by the parameter M LyC ffT2l . The earthquake times are denoted by the random 
variable £j, where the index i counts the natural time (i.e. the event order), and the interevent times 
are denoted by r< = U+i — U. Finally, A will denote the estimate of A from data, where A stands 
for a function or a distribution parameter. We also use the following abbreviations: PDF for the 
probability density function, CDF for the cumulative distribution function, ITD for the interevent 
times distribution and CSS for the Cretan seismic sequence. 

B. Earthquake interevent times 

The distributions of the magnitudes and return intervals of extreme events are research topics 
that attract significant interest [fT3l . Extreme events are usually distributed in both space and time, 
but herein we assume that the spatial dependence can be ignored to simplify the analysis. Let us 
further assume that an extreme event corresponds to excursions of a random function X(t), where 
t is the time index, to values above a specified threshold z q . The classical extreme value theory 
focuses on the probability distributions of such extreme events lfT4l . The Fisher- Tippet- Gnedenko 
theorem |[T5l [T6l shows that the distribution of M n : = max(X 1 , . . . ,X n ), where the X^i = 
1, ... ,n are independent and identically distributed (i.i.d.) variables, is given by Generalized 
Extreme Value (GEV) distributions. The GEV include the Gumbel, reverse Weibull and Frechet 
distributions. The type of GEV obtained depends on the tail behavior of the probability distribution 
of the Xi. Similar distributions, but with reversed supports, are obtained for minima. Extension of 
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the classical CLTs (which involve deterministic scalings of the random variables), to randomized 
CLTs that enable stochastic scaling transformations have been recently proposed in ifTTll . 

If the earthquakes in a fault system occurred randomly, the event times should be uniformly 
distributed. Uniform distribution over an unbounded time domain implies Poisson statistics as lu- 
cidly explained in IfTTll . The Poisson model leads to an exponential distribution of interevent times. 
It has been proposed that the earthquake times follow the Poisson distribution if foreshocks and 
aftershocks are removed lfT8l . Different declustering approaches have been proposed to isolate 
main events from aftershocks and foreshocks. Nevertheless, these approaches are not based on 
fundamental principles, and the results of their application on seismic data vary widely. In addi- 
tion, a recent rigorous statistical analysis casts doubt on the validity of the Poisson model even for 
declustered data [6|. The periodic model supports that the intervals between characteristic earth- 
quakes are approximately constant lfT9l . This is in contrast with the observations collected from 
large faults globally lr5ll20ll2D. 

Spatial and temporal inter-dependence of seismic events can explain deviations from the Pois- 
son and the periodic models. Several research papers propose that earthquakes are self-organized 
systems, perhaps near a critical point flTJIEllTTJ or systems near a spinodal critical point ll4l l22l. 
Both cases imply the emergence of power laws in the system. Bak et al. [7] introduced a global 
scaling law that relates earthquake interevent intervals with magnitude and distance between earth- 
quakes. These authors analyzed seismic catalogue data from an extended area in California that 
includes several faults over a period of 16 years (ca. 3.35 x 10 5 events). They observed power-law 
dependence over 8 orders of magnitude, indicating correlations over a wide range of interevent 
intervals, distances and magnitudes. Corral and coworkers |j8] 4T0ll23l introduced a local modifica- 
tion of the scaling form so that the interevent time probability density function (PDF) follows the 
universal expression f T (r) ~ A/ (At), where /(r) is a scaling function and the typical interevent 
time f is specific to the region of interest. 

Saichev and Sornette IfTTl [24) generalized the scaling function incorporating parameters with 
locally varying values. Their analysis was based on the mean-field approximation of the in- 
terevent times PDF in the epidemic-type aftershock sequence (ETAS) model. ETAS incorporates 
the Gutenberg-Richter dependence of frequency on magnitude, the Omori - Utsu law for the rate 
of aftershocks, and a similarity assumption that does not distinguish between foreshocks, main 
events and aftershocks (any event can be considered as a trigger for subsequent events). 

In a different research vein, several studies of earthquake catalogues and simulations show that 
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the Weibull distribution provides a better match of empirical interevent times distributions than 
the Poisson model B51 125U3TI . The arguments supporting the Weibull distribution are based on 
empirical studies and extreme value theory. In particular, since the interevent times are determined 
by minima of the shear strength in the Earth's crust, the standard Weibull model is a good candidate 
for their distribution. In contrast, the Gumbel distribution for minima has negative support and the 
Frechet distribution has an unbounded support. 



C. Aims and outline of this paper 

In this paper we propose a stochastic stick - slip model that links the shear strength (in the 
following referred to as strength for brevity) distribution of faults in the Earth's crust, the stress 
accumulation - relaxation process in the crust due to tectonic motion, and the earthquakes in- 
terevent times distribution. The model is formulated at the single-fault scale and is then extended 
to a system of faults by constructing a composite strength distribution. 

A prototype physical model of stick - slip motion along single faults is the Burridge-Knopoff 
(BK) model Il32l - |34l . This model consists of a system of coupled differential equations represent- 
ing the motion of n point masses linked with elastic springs and subject to a velocity-dependent 
friction force, which is responsible for the slip instability. In contrast to the BK model, the stochas- 
tic stick - slip model proposed herein is phenomenological, since the time between events is deter- 
mined by a heuristic stress accumulation function. We assume that the main stochastic component 
is due to the variations of the fault shear strength (or the strength across different faults in a sys- 
tem). Stochastic aspects of the stress accumulation function are not explicitly investigated in the 
following, but they obviously deserve further research. 

The remainder of the paper is organized as follows: In section|II]we review the Weibull distri- 
bution and its applications to earthquake interevent times. In Section [III] we propose a stochastic 
stick - slip model for single faults and show that it admits interevent times distributions, includ- 
ing the Weibull, that depend on the time evolution of stress accumulation. In particular, we show 
that the Weibull is an admissible interevent times model, if (i) the crust strength distribution is 
Weibull with stationary and homogeneous parameters, and (ii) the stress accumulation increases 
with time as a power-law. We also show that deviations from the Weibull can result due to spatial 
non-homogeneity of the crustal strength parameters and by imposing finite magnitude thresholds 



on the seismic sequence. In section IV we generalize the model to interevent times for a system 
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of faults. In Section [V] we investigate the ITD for a microearthquake sequence from the island 



of Crete (Greece) in relation to the proposed stick - slip model. Finally, Section VI involves a 
discussion, conclusions, and topics for further research. 

H. THE WEIBULL DISTRIBUTION AND EARTHQUAKE INTEREVENT TIMES 

A. Properties of the Weibull distribution 

The Weibull CDF, F T (r), determining the probability that the time between two consecutive 
events is less than or equal to r > is given by 



F t {t) = 1 — exp 



, m-i 

T 



(1) 



where t s is the scale parameter and m > is the Weibull modulus or shape parameter. The PDF 
is defined by / T (r) = dF T (r)/dr and is given by the equation 

/r(T) = ™ fr)™" 1 e -(-r. ( 2) 

The survival function is the complementary cumulative probability function, i.e., 

R(t) = 1-F t (t). (3) 

In the Weibull case, R(r) is the stretched exponential function R(r) = e . The function 
R(t) represents the probability that no seismic event has occurred within time r since the last 
event. Shape parameter values m < 1 lead to a diverging density at r = 0, and an almost 
exponential decay of / T (r) as r increases. For m > 1 the PDF develops a single peak that 
becomes sharper with increasing m. 

The hazard rate or failure rate is the rate of change for the probability of a seismic event, if 
time t has elapsed since the last event. It is expressed by 



/t(t) m ( r 



m—l 



H(r) = - ^ = - " • (4) 
\-F t [t) t s \t s J 

Data are graphically tested for Weibull dependence using the Weibull plot. The latter employs 
the fact that the double logarithm of the inverse survival function, 

$(r) -loglogiT^r), (5) 

satisfies the straight line equation, $(r) — m log(r) — m log(r s ) with slope equal to the Weibull 
modulus. If the data are drawn from the Weibull distribution, $(r) obtained from the empirical 
CDF, F t (t), is approximately a straight line. 
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B. The Weibull interevent times distribution 

The Weibull distribution was investigated in 11251 - 1271 for large earthquakes at six subduction 
zones over the globe. In Il28l the ITD of a sequence of 12 paleoearthquakes on the San Andreas 
fault was investigated. Yakovlev et al. [|29l simulated million-year-long catalogues of earthquakes 
on major strike-slip faults in California. They found that the Weibull distribution fits the interevent 
times better than the lognormal and inverse Gaussian distributions. 

Abaimov et al. J51 |35j concluded that the Weibull is a good model for real and simulated 
large-magnitude earthquakes on the San Andreas fault, and for a microearthquake sequence at a 
nearby site. These authors also emphasize the behavior of the hazard rate function of the Weibull 
distribution |[36] - [38l . At least for large-magnitude earthquakes, H{r) is expected to increase with 
the interevent time. The exponential distribution has a constant H[t) indicating lack of memory 
between events, while for the lognormal distribution H(r) decreases with the interevent time. The 
Brownian-passage time distribution H(r) tends to a constant with increasing interevent time. Of 
the various distributions considered as models of earthquake interevent times, an increasing H(r) 
with time is exhibited only by the Weibull and gamma distributions (if m > 1). In addition, 
Abaimov et al. find evidence for the Weibull distribution in numerical solutions of the slider- 
block model introduced by Burridge and Knopoff 113211391 . 

Robinson et al. [40] simulated approximately 500 000 earthquakes of magnitudes between 
3.8Mw and 6.6Mw (moment magnitude scale), over a period of two million years for faults in 
the Taupo Rift in New Zealand. These authors employed a synthetic seismicity model that is 
based on the Coulomb failure criterion and uses empirical data pertaining to the number of faults, 
fault lengths, and long-term slip rates. They found that a three-parameter Weibull distribution fits 
the interevent times of large earthquakes on most faults. The three-parameter Weibull survival 
function is given by R{r) = e~( r ~ T °)" l / r ™ ', where r is the location parameter. A finite r implies 
that the PDF vanish for r < r . According to the analysis in 11401 , earthquakes within a normal 
fault system are correlated on a rift- wide scale over time periods of w 3kyrs. 

Santhanam and Gantz [HTTl investigated the return intervals of a random function X(t), i.e., 
the times between consecutive excursions of X(t) above a given threshold. They found that if 
X(t) has long-range memory (i.e., power-law decay of correlations), the return intervals follow 
the Weibull distribution. Their mathematical formulation provides support for the Weibull model. 
Nevertheless, their approach is not based on the standard laws of seismology (i.e., Gutenberg- 
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Richter and Omori's law), and the values of m admitted are < m < 1. In contrast, analysis 
of interevent times from seismic catalogues and simulations demonstrates that values m > 1 also 
occur. In particular, for large earthquakes it is believed that the hazard rate function increases with 
the time since the last event, implying m > 1. 



TABLE I: Weibull modulus (m) values for earthquake ITD from published research based on 
empirical and synthetic data from various tectonic settings and system sizes. The magnitudes in 
column 4 are based on the moment magnitude scale (Mw) [|42l|. 



# of events 


m 


Time span 


Magnitude (Mw) 


Type of data 


System "Size" 


Location 


Reference 


7 


3,21 


147 yrs 


> 6 


Real 


Single Fault 


r~i 

SArM Pallett Creek 


® 


13 


2,01 


ft2000 yrs 


> 7 


Real 


Single Fault 


SAF/Wrightwood 


ED 


13 


2,91 


ft: 5 years 


r~| 

Mean=1.3fP 


Real 


Single Fault 


SAF/San Juan Batista 





4606 


1,97 


r-i 

i mJ c 


> 7.5 


Simulation 


Fault system 


North SAF/S.S. faults 


ED 


5093 


1,87 


1 Ma 


> 7.5 


Simulation 


Fault system 


South SAF 


ED 


2612 


1,71 


1 Ma 


> 7 


Simulation 


Fault system 


Hayward / California 


(29) 


8174 


1,42 


1 Ma 


> 6.8 


Simulation 


Fault system 


Calaveras / California 


(29) 


1913 


1,32 


1 Ma 


> 6.7 


Simulation 


Fault system 


San Gabriel / California 


|29| 


1075 


1,7 


1 Ma 


> 7.2 


Simulation 


Fault system 


Calaveras / California 


(29) 


7 


2,9 


147 yrs 


> 6 


Real 


Fault system 


SAF / Parkfield 


(45) 


7 


1,5 


« 2000 yrs 


> 7 


Real 


Fault system 


SAF/ Pallett Creek 


(U 


510 s 


0.82-10^] 


2 Ma 


3.3-6.8 


Simulation 


Fault system 


Taupo Rift / New Zealand 


(40) 


12024 


0,91 


82 months 


> 4 


Real 


Arej=] 


Okinawa/Japan 


ED 


13678 


0,79 


82 months 


> 4 


Real 


Area 


Chuetsu/Japan 


ED 


12024 


1,09 


82 months 


> 3.5 


Real 


Area 


Okinawa/Japan 


ED 


13678 


0,85 


82 months 


> 3.5 


Real 


Area 


Chuetsu/Japan 


ED 


12024 


1,43 


82 months 


> 3 


Real 


Area 


Okinawa/Japan 


ED 


13678 


1,08 


82 months 


> 3 


Real 


Area 


Chuetsu/Japan 


ED 


12024 


1,75 


82 months 


> 2 


Real 


Area 


Okinawa/Japan 


ED 


13678 


1,77 


82 months 


> 2 


Real 


Area 


Chuetsu/Japan 


ED 




1,7 


70 years 


6.9-8.4 


Real 


Area 


Japan 


ED 


8 


2,3 


1000 yrs 


7.9-8.4 


Real 


Subduction margin 


Nankai Trench / Japan 


ED 


16 


2,9 


300 yrs 


7.8-8.4 


Real 


Subduction margin 


Hokkaido-Kurille Trench / Japan 


ED 


11 


9,5 


100 yrs 


7.8-8.7 


Real 


Subduction margin 


Aluetian Trench / Alaska 


ED 



a SAF: San Andreas Fault, California. 
b Analysis of data from |44l . 
c lMa=l Million years 

d In ll40l the data are fitted to the three-parameter Weibull distribution. 

e This refers to large areas that may include several fault systems and subduction margins. 



Table previews published estimates of the Weibull parameters for the interevent times of various 
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earthquake sequences, along with information pertaining to the magnitude, the size of the fault 
system, the location, and the duration of the seismic sequence. Note that values of m > 1 prevail. 
In addition, estimates of m obtained from records containing a small number of events tend to be 
higher than those from large sample sizes that lie in the range 0.79 — 1.97. We believe that this 



tendency is partly due to the impact of the magnitude threshold on the ITD (also see section III B 



below). Also note that the data cover very different system sizes, from single faults to areas 
containing many fault systems. 



III. A STOCHASTIC STICK - SLIP MODEL FOR SINGLE FAULTS 

An overview of the physical processes associated with earthquakes is given by Kanamori and 
Brodsky [|46l|47]|. We develop our formalism based on this conceptual framework. Specifically, 
we establish connections between the shear strength distribution of the Earth's crust, the dynamic 
process driving earthquake generation, and the distribution of interevent times. Tectonic earth- 
quakes occur on fault planes due to a dynamic stick-slip process that locally accumulates stress 
caused by the plate motion. This stress accumulation eventually leads to crust failure and stress 
relaxation when the local crustal strength is exceeded. The phase of stress accumulation (loading 
phase) is followed by a phase of rapid stress relaxation, cf. Fig. 3(b) in Il47ll . which corresponds 
to the slip events. The stress accumulation - relaxation process is cyclically repeated. This model 
has been applied with constant stress accumulation rate and constant maximum, er max , or residual, 
cr res , stress to explain the recurrence of large earthquakes. If both the maximum and residual stress 
are constant, the model predicts periodic behavior; if the initial stress is constant, the model pre- 
dicts the time of the next event; if the final stress is constant, the model predicts that the longer the 
interevent time, the larger the magnitude of the following event ||48i 



A. On the distribution of crustal strength 

Earthquakes are typically localized on faults in the Earth's crust; hence, the strength of these 
structures and the applied tectonic loading determine the interevent times on single faults. Since 
the crust is composed of brittle material (rock), its strength is expected to follow the Weibull prob- 
ability distribution [|49ll . The Weibull distribution was derived in the framework of weakest-link 
theory, founded in the studies of Gumbel and Weibull on extreme- value statistics 11501 . This theory 
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addresses the strength of brittle and quasibrittle materials in terms of the strengths of "represen- 
tative volume elements" (RVEs) or links [|5T1455l . The material fails if the RVE with the lowest 
strength breaks, hence the term "weakest-link". The concept of links is straightforward in simple 
systems, e.g., one-dimensional chains and fibers. In the case of more complicated structures, the 
links represent critical units that control the failure of the entire system. 

Experimental studies on the strength of geological materials, such as rock (granite) under var- 
ious types of loading (compressive, bending) provide evidence for the validity of the Weibull dis- 
tribution It56l 15711 . At the same time, experimental measurements show that the crustal strength has 
a systematic dependence on the depth h ll58T[60ll . Hence, for a single fault the following Weibull 
CDF is a useful approximation of the strength distribution at depth h: 

F s [S;S 8 (h)] = l-e-{*™> , (6) 

where S s (h) is the strength scale and m s is the strength Weibull modulus. Typical values of m s 
for laboratory measurements of rock samples range between 3 and 30 It56ll57l . 

The seismic events on a fault are distributed across different depths in the crust. The systematic 
dependence of crustal strength measurements on depth ll59l - l6H implies that S s (h) increases con- 
tinuously with h. This agrees qualitatively with the decreasing number of seismic events registered 
with increasing depth (c.f. Section |V| below). 

The depth- averaged effective fault strength distribution is then given by 

1 f h2 „ -l^- 



F s(S) = 1 - j^—f^ / dhe-^) , (7) 



In 



where hi is the minimum and h 2 the maximum focal depth of the earthquake events. In Ap- 
pendix [A] we evaluate the effective distribution for a linear dependence S s (h) = S + qh, where 
S s (hi) = S and S s (h 2 ) = S + q h 2 . It is shown that 

FS(s) = 1 - — Kr f 1 duu- {l+1/ms) e~ bu , b = s m % (8) 



where s = S/S is the normalized strength, 5h = 2(h 2 — hi), X s = q5h/S is the strength 
variability coefficient, and iii )2 = (i^br) ■ ^ n Appendix |a| we also derive an explicit expression 
for the integral valid for m s > 1 and a convergent infinite series that is valid for m s < 1. 

In Fig.[T]we show the numerically integrated CDF Fg (s), as well its difference from the Weibull 
CDF,F (s) = 1 — e~ sms , i.e., AFg (s) = F§(s) — F (s). Larger m s values lead to more pronounced 
AFg (s) . Weibull plots of the effective strength distribution, shown in Fig. [5[ reveal that for m s < 1 
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FIG. 1: F*(s) (top left) and AF*(s) (top right) for X s = 0.7 and < m s < 1. F*(s) (bottom 
left) and AF*(s) (bottom right) for \ s = 0.7 and 10 > m s > 1. 

the linear dependence of $*(s) = log log 1 _p-^ s - ) is maintained after the averaging, while $*(s) 
deviates from the linear (Weibull) dependence for m s > 1 . 

B. Relation between crustal strength and ITD 

The stochastic stick - slip model we propose incorporates non-uniformity of the crustal strength 
and links the stress accumulation and relaxation processes with the strength distribution and the 
earthquake interevent times. 

Let us assume that tj G [0, tj] measures the time during the i-th loading phase, and that the stress 
increase is determined by a(ti) = 4>(U), where the loading function 0(t) is an increasing function 
of time. In general, the parameters of 4>(t) vary randomly between different phases. For simplicity, 
we first assume that <j>(t) is a deterministic function, the relaxation time is negligible, and the stress 
relaxes to a zero residual value after a seismic event. Let the random variable a(ti) correspond 
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log(s) log(s) 

(a) $(s) for < m s < 1. (b) $(s) for 10 > m s > 1. 

FIG. 2: Weibull plots of the effective fault strength distribution based on Eq. @. 

to the stress at failure for the seismic event E{. Hence, a maXj j := a(ti) is equal to the strength, 
Si, of the crust at the particular location and time, and we can assume that it follows the Weibull 
distribution. The time of failure is given by tj = 0~ 1 (S' i ), where represents the inverse of 

4>(t). Once the event takes place, the stress is relaxed and the next event E i+ i occurs when 
Si+i = <j(i i+ i). The interevent time between events E { and E i+ x is given by t i+ i = 0~ 1 (S'j + i). In 
general, if the crustal shear strength S is viewed as a random variable, it is related to the interevent 
times random variable, r, by means of 5* = 0(r). 

If fs(S) is the PDF of the crust strength, and 0(r) is a differentiable, monotone function, the 
corresponding PDF of the interevent times is obtained by means of Jacobi's theorem for univariate 
variable transformations as follows: 

/ r (r) = / s (0(r))|0 / (r)|, 

where 4>'{r) is the derivative of the loading function. In particular, if the crustal strength follows 
the Weibull distribution with scale parameter S s and modulus m s , the ITD has a PDF that is 
determined by the equation 
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C. Stress accumulation scenarios 



1. Linear stress accumulation 



a. Zero residual stress: If = a r, where a is the stress accumulation rate, it follows 
from Eq. © that the ITD is the Weibull distribution with the CDF of Eq. Q with r s = S s /a and 
Weibull modulus m = m s . In this scenario the stress relaxes to zero after each seismic event, c.f. 



Fig |3aJ Since the stress accumulation rate is independent of the strength, t s and m may also vary 
independently. 

b. Finite residual stress: If the stress relaxation process terminates at a non-zero residual 



stress, cr rcs , the loading function is given by = a res + at. This scenario, illustrated in Fig. 3b 
is equivalent to the elastic rebound theory developed by Reed [|2l l62l| : if a fault with a shear 
modulus G and width 2b is sheared at constant tectonic velocity u , then a = Gu /2b. Since 
the residual stress does not cause failure of the crust, the crustal strength is expected to have a 
minimum value S min > <r res . Hence, the strength distribution corresponds to a three-parameter 
Weibull, the CDF of which is given by Eq. ([T]) with a replaced by a — S min . In this case, based on 
Eq. ® the ITD becomes 

fAr) = n ( L^Y"" 1 e -F=M- (10) 

where t\ oc = (S min — a ICS )/a is the time location parameter. If S min = cr res the two-parameter 
Weibull model is recovered. In addition, if S m i n — a TCS « S s , the two-parameter Weibull is an 
accurate approximation. A value of Ti oc that is not negligible compared to r s supports the use of 
the three-parameter Weibull ITD model in [|40l . 

c. Finite relaxation time: The analysis above is not substantially modified if the relaxation 



has a finite duration t* +1 (as shown schematically in Fig. 3c), and is governed by a decreasing 
function (j) re \(t), since t* +1 = (0 r ei) _1 (S i+ i) is also determined from the strength S i+ i. For exam- 
ple, assuming linear loading and relaxation relations, it follows that the stress at failure is given by 
o~(ti+i) = ott i+ i and a(t i+1 ) = a rel f* +1 . Hence, the total interevent time (measured between the 
zero stress values of consecutive phases) is given by t i+ \ + t* +1 , which is equivalent to replacing 
a^ 1 with oT x + a~ c \. Hence, in the following we simply use a without loss of generality. In 
the case of aftershocks triggered by a main shock, the relaxation may involve a non-monotonic 
evolution of the stress toward its residual value. Then, a rc i should be considered as an effective 
relaxation rate and modeled as a random variable. 
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FIG. 3: Schematic drawings of stress accumulation scenarios, (a) Linear loading function with 
zero residual stress and instantaneous relaxation, (b) Linear loading function with finite residual 
stress (c) Linear loading function with finite relaxation time, (d) Nonlinear loading function, (e) 
Linear loading function with stochastic rate, (f) Stochastic loading function. 

2. Nonlinear stress accumulation 



The following scenarios assume that the relaxation time is negligible and that the residual stress 
is zero. 

a. Logarithmic stress accumulation: This case is governed by a loading function 0(r) = 
a* logfY/Vd), which implies sublinear stress evolution. According to Eq. (|9]), logarithmic loading 
leads to the log-Weibull ITD 



T 



log 



m s — l 



r* log(T/r d ) V 



(ID 



An equivalent expression has been investigated in ||3~0ll3"TL where the authors employ a linear 
mixture of the Weibull and log-Weibull distributions to model the ITD for earthquakes in different 
tectonic settings. As we have shown above, the log-Weibull is justified for a sub-linear stress 
accumulation. 



b. Power-law stress accumulation: This scenario, depicted in Fig. 3d is governed by the 
loading function 4>(i) = ar 13 , which leads to the Weibull ITD 
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m— 1 

r ' 



f T (r) = m - e-K-s) , (12) 



where m = m s (3 and r s = {Ss/a) 1 ^ . In the following, we do not distinguish between the linear 
and the power-law stress accumulations models since they both lead to a Weibull ITD. 

Note that the relation m = m s (3 enables an indirect estimation of the loading exponent (3 
from the ITD Weibull modulus m, obtained from the analysis of the seismic event times, provided 
that the strength Weibull modulus is estimated from laboratory measurements. For example, if 
m s = 10 and m = 0.7, it follows that (3 = 0.07. Let us also use the transformation a = A/tq , 
where A is a constant with units of stress and t is an arbitrary time constant. Then, it is possible 
to empirically determine A. For example, if S s is known, e.g., S s = 250 MPa (mega Pascal), 
r s = lday while r is arbitrarily set to r = lsec, it follows that A = S s / (r s /r ) /3 ~ 87.57 Pa. 

c. Stochastic stress accumulation: In general, the statistical properties of the crust change 
over geological time fl47). Such changes on individual faults may be due to progressive damage 
caused by ongoing seismic activity Q. Hence, the parameters of the strength distribution may 
exhibit variations over time. In addition, the stress accumulation process may exhibit stochastic 
behavior over different time scales. For example, in the linear stress accumulation scenario the 



rate may fluctuate, as shown schematically in Fig. 3e More generally, the stress accumulation 



may have a complex dependence with intra-phase variations of the accumulation rate as shown in 



Fig. 3f In such cases, the fault is characterized by the effective ITD: 



F;(t;6*)= f U d6p(0)F T (T;8), 

JBd 

where = (9\, . . . , 6k) is a parameter vector with joint PDF p{0) that includes both strength 
and stress accumulation parameters. The function p{6) is non-negative and integrable. If all the 
components of 6 have a finite support and F T (r; 6) is a continuous function of 0, it is possible to 
iteratively apply the first mean-value theorem of integration, which leads to 

F;(t; 0*) = F T (r; 6\{r), 0*(r, &{),... 6* k (r, &U)) (13) 



If the dependence of the effective parameters 6* on r is weak, the effective CDF is of the same 
functional form as F t (t; 6), with 6 replaced by the effective parameter vector 0* whose compo- 
nents are inter-dependent. 
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D. Magnitude dependence of ITD on single faults 

In the following, we refer to the interevent times for the complete (i.e., including all magni- 
tudes) seismic sequence as primitive interevent times. The main assumption in the analysis is that 
the primitive ITD conforms to the Weibull distribution. 

1. Impact of magnitude threshold on sampled strength distribution 

To link the model of repeated stress accumulation and relaxation phases with the ITD, we 
have assumed that all the seismic events are included, regardless of their magnitude. Then, the in- 
terevent times sample the entire crust strength distribution. In practice (see Section[V]), one focuses 
on seismic events over a threshold magnitude that varies depending on the area, the resolution of 
the instrumental network, and the goals of seismic risk assessment. Often the threshold magnitude 
is chosen as the magnitude of completeness above which all events are resolved by the observation 
system [63]. 

The stress drop caused by an earthquake is related to the earthquake magnitude by means of 
empirical, monotonic relations ll64ll . Then, assuming that a Tes is constant, the threshold magnitude, 
Ml, c , corresponds to a unique crustal strength value S c . Hence, the proportion of earthquakes 
exceeding in magnitude M Lc is 1 — Fg(S c ). The sampled PDF for events with S > S c is given by 
the following equation 

fU(S) = ! _ /s ^ c) W-S e ), (14) 

where $(•) is the unit step function and the denominator normalizes the PDF. The removal of low 
strength values implies that the respective $(S) has a concave lower tail, even if the effective fault 
distribution Fg (S) is Weibull. In this respect, the impact of coarse-graining by means of a finite 
threshold is similar to that of non-resolved, low-magnitude events. 

In Fig. [4] we show the impact of coarse-graining on a sample of Weibull random numbers with 
S s = 100 MPa (which is a typical average value of normal-fault strength [58]) and m s = 0.7 with 
cutoffs S c = 300 KPa and 3 MPa. In spite of the curvature of $(S f ) below the threshold, the linear 
part of both ®(S) plots has the same slope, i.e., the same Weibull modulus. Nevertheless, the 
estimated modulus ll65ll based on the truncated sequence increases with S c ; namely, m s « 0.74 
for S c = 0.3 MPa and m s ps 0.83 for S c = 3 MPa. Hence, an increasing S c leads to an increase of 
the estimated Weibull strength modulus. 
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FIG. 4: Plots of &(S) derived from a sample of 10.000 Weibull random numbers from a 
distribution with S s = 100 MPa and m s = 0.7 using two different thresholds (S c = 300 KPa, 3 

MPa). 



The times between events exceeding Ml jC extend to several stress accumulation and relax- 
ation cycles. We assume a stress accumulation phase between events with Ml > Ml >c , such that 
stress increases non-monotonically between two events from a low residual value to the strength 
threshold that corresponds to M L . Then, we can replace the non-monotonic dependence with 
an "effective" monotonically increasing stress accumulation function; e.g., in the linear case, a 
is replaced by a e g. Consequently, the ITD for finite M L jC is also determined from Eq. ( p"4] ). In 
general, a e fr should be treated as a random variable that fluctuates depending on the number of 
sub-threshold events contained between the two supra-threshold events. A power-law stress ac- 
cumulation function transfers to the ITD the lower-tail curvature of the crustal strength resulting 
from the truncation of above-threshold events. 

The interevent times between events exceeding Ml )C can be treated as sums of a randomly 
varying number of primitive times. To our knowledge, very little is known about sums of Weibull 
random numbers: a review of known results for sums of random variables reports the absence of 
such results for the Weibull distribution [66J. Recently, convergent series expansions for the PDF 



2. Impact of magnitude threshold on ITD 
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that governs the sum of Weibull random numbers were obtained in ll67l . The situation is more 
complicated for earthquake interevent times, since the number of summands, which is determined 
by the excursions of the dynamic stress above S c , fluctuates. This coarse-graining operation com- 
pletely eliminates strength values below S c , in contrast with the summation of a fixed number of 
variables that shifts the mode of the distribution but does not eliminate the weight of the PDF near 
zero. 



IV. ITD OF FAULT SYSTEMS 

Above, we focused on the one-body ITD problem by concentrating on single faults. Neverthe- 
less, the seismic behavior of a given area is a many -body problem that involves multiple faults. 
The interevent times of a fault system are not directly obtained from the interevent times of the 
individual faults in the system. Let us assume that a fault Fa hosts two seismic events at times 
t Fa and tj^ leading to an interevent time r Fa = tj+j — t^ a . In addition, let fault Fb host two 
seismic events at times tf b and t[K, with an interevent time rf b = tfh — tf b . Next, let us assume 
that tf a < tf b < tf b x < t?+ x . Then, the respective interevent times for the two-fault system are 
tf b - tf a , tf b x - tf b , tf^ - tf b , which can not be obtained from r[ a and rf b . 

Nevertheless, we can extend the single-fault interevent-times expressions to fault systems by 
assuming that all the faults are subject to a uniform stress accumulation process. Then, the same 
loading function applies to the entire system, and the ITD is determined from the composite 
strength distribution of the system ll68l . Assuming that the system involves N faults characterized 
by K different effective strength distributions Fg(S;6i), % — 1, . . . , K and that Mj is the number 
of faults that follow the probability distribution F£(S\ Oi) (i.e., M x + M 2 + . . . M K = N), the 
composite strength distribution of the fault system is given by the following superposition 

^0S) = Z>*s Pi = ^,i = l,..-K. (15) 

i=i 

A homogeneous fault system comprises faults that share the same strength distribution, Fg (S; 6). 
Then, the composite strength distribution is given by Jg(S) = F^(S; 6). 
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A. System with bimodal strength distribution 

Let us consider a system that involves faults governed mainly by two Weibull strength distribu- 
tions with different parameters, so that K\ faults follow the first distribution, while K 2 = N — K± 
faults follow the second. The composite strength distribution of the system is a bimodal Weibull. 
Let us assume that the survival function is given by the bimodal expression 

f S V m l ( S \ m 2 

1 - jr*(5) =pe UJ +(l-p)e W , (16) 

where p = K\/M. This mixture model can lead to the "M-smile" histogram sometimes observed 
in the analysis of earthquake interevent times, e.g. Il69~ll . 

In Fig. [5] we show the histogram and corresponding empirical $(r) for a set of 3000 interevent 
times generated from the bimodal Weibull distribution with p = 2/3, m 1 = 0.7, m 2 = 1.8, 
T\ = 1.5 10 3 (sec), and r 2 = 10 5 (sec). The signature of bimodality is apparent as a dip in the 
histogram and as a saddle point in the Weibull plot. The parameters of a bimodal ITD can be 
estimated using the expectation-maximization algorithm lf70l . 

The concept of a bimodal ITD has been proposed in [69], where it is invoked to represent the 
mixing of correlated events (aftershocks) and uncorrelated (main) events. In our model, there is 
no a priori distinction between aftershocks and main events. However, it is conceivable that main 
events lead to temporary changes of the stress accumulation or the crustal strength, which imply a 
distinct component in the ITD linked with aftershock activity. 

B. Magnitude dependence of ITD for fault systems 

For a single fault we showed above that the lower tail of the sampled strength PDF is cut-off 
(c.f. Fig [4}. We expect that the abrupt change in the slope of $(S) at S ps S c is reduced in 
seismic data from fault systems due to non-homogeneity of the strength parameters and possible 
non-stationarity caused by fluctuations in the Weibull parameters over time. 

In the following, we assume that the fault-system ITD for events exceeding M L:C is approxi- 
mated by the Weibull. The total duration of the complete seismic sequence is T = J2f=i T i- Let 
t\ and T2 denote sample average times corresponding to truncated seismic sequences obtained for 
Mi l < Ml ; 2, respectively. We assume 1\ and T2 to be accurate estimates of the ensemble means 
of each sequence, i.e., v[/t^ w E[ti]/E[t2]. If N t is the number of events with magnitude above 
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Interevent times (sec) 



FIG. 5: Frequency histogram (log-log plot) and Weibull plot (inset) for 3000 random numbers 
simulated from a bimodal Weibull mixture model (see text for parameters). 



M Lti , then T { « T/N { , We assume that the Gutenberg - Richter scaling, N 1 /N 2 = lO" 6 ^. 1 "^. 2 ), 
is valid for the system. Then, it follows that 

?t E[n] 

where p M = b log 10. Since for seismically active regions b e [0.5 — 1.5], it follows that pu G 
[1.15-3.45]. 

The mean value of the Weibull distribution is given by E[r] = r s r(l + l/m). Hence, r Si2 /T s ,i = 
E[r 2 ] T(l + l/mi)/E[n] T(l + l/m 2 ). Then, based on Eq. ([17) it follows that 

Ts,2 _ ^1 + 1/mO „- PM (M L , 1 -M i , 2 ) 



r,,i T(l + l/m 2 ) 



e ~PM Ml i2 ) (IS) 



Equation ( fT8| ) also holds for single faults that obey Gutenberg - Richter scaling and have a Weibull 
strength distribution. 
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V. ANALYSIS OF CRETAN MICRO-EARTHQUAKE SEQUENCES 

The seismic data in the Cretan seismic sequence (CSS) investigated below are from Becker et 
al. irTDl . The CSS resulted from tectonic activity generated at the Hellenic subduction margin, 
where the African plate is being subducted beneath the Eurasian plate; this is the seismically 
most active region in Europe. More than 2 500 local and regional micro-earthquake events with 
magnitudes up to 4.5 Ml (Richter local magnitude scale) occurred during the time period between 
July 2003 and June 2004. The micro-earthquakes were accurately recorded by an amphibian 
seismic network zone onshore and offshore Crete. The network configuration consisted of up 
to eight ocean bottom seismometers as well as five temporary short-period and six permanent 
broadband stations on Crete and smaller surrounding islands (e.g. Gavdos). The magnitude of 
completeness varies between 1.5 Ml on Crete and adjacent areas and 2.5 Ml at around 100 km 
south of Crete. Most of the seismic activity is located offshore of central and eastern Crete (see 
Fig. [6]). The repeat times between successive earthquake events range from 1 sec to 19.5 days. 



A. Exploratory analysis 

The exploratory analysis of CSS includes all the events in the catalogue (magnitudes > IMl). 
There is no discernible main shock in the sequence, hence declustering of the data is not consid- 
ered. We use the concept of natural time to measure the ordering of the event times [[7211731 : if t k 
is the time of the fc-th event, the interevent time series in natural time is given by r k := t k+1 — t k . 



The CSS interevent times series is shown in Fig. 7a , which exhibits isolated large peaks separating 



clusters of almost continuous seismic activity. The distribution of the earthquake focal depths is 
shown in Fig.|7bJ Assuming a uniform tectonic stress over depth, the observed declining trend [|74l 
agrees qualitatively with the reported increase of the crust strength with depth ll58l 



The histogram of the magnitudes is shown in Fig. 8a The Gutenberg - Richter law is expressed 



as log 10 N(M L > M) = a — bM, where iV is the number of events with magnitude greater or 



equal to M. This implies an exponential decrease of N with the magnitude, while Fig. 8a shows a 
non-monotonic dependence of iV on M. This is due to the resolution problem, namely the inability 
to observe all events with magnitudes below a critical threshold that is typically around 2 Ml and 
m 2.2Ml for the CSS. Setting the critical threshold at 2.2 Ml, the histogram of event frequency 



versus magnitude is shown in Fig. 8b The maximum likelihood estimate of the Gutenberg 
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Richter exponent is b w 0.85, i.e., within the typical range [0.5 — 1.5] for seismically active 
regions. 

The Weibull plot for CSS is shown in Fig.[9j The slope of the best-fit (minimum least squares) 



straight line is rh = 0.66. Note that based on the analysis in III B , since m < 1 we do not 
expect significant modification of the Weibull due to the depth dependence of S s . Nevertheless, 
the lower tail of $(r) is lighter than the Weibull while its upper tail is heavier. The lower tail 
begins roughly at $(ti) ~ —2.5 and the upper tail near $(r 2 ) ~ 1.4. By inverting Eq. (|5]), we 
determine that F t (ti) ~ 0.08 and F t (t 2 ) ~ 0.98. The lower-tail behavior, which involves about 
8% of the data, can be attributed to unresolved events (i.e., with magnitudes below the magnitude 
of completeness). The upper-tail behavior, which involves around 2% of the data, can be partly 
explained due to the same effect, since failure to observe certain seismic events leads inadvertently 
to interevent times that are higher than in reality. The discrepancy could also be due to genuine 
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FIG. 7: (a) Sequence of interevent times: the horizontal axis counts the order (i — 1, . . . , N) of 
the interval Tj = ti + \ — ti, and the vertical axis measures Tj. (b) Histogram of the focal depths of 

earthquake events. 
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FIG. 8: (a) Frequency histogram of events (complete CSS) versus magnitude, (b) Histogram of 
event number versus threshold magnitude and best-fit line to Gutenberg - Richter model. 



departure of the upper ITD tail from the Weibull behavior, since the CSS data set involves a system 
of faults with a composite crustal strength distribution, c.f. Eq. ([15]). 



25 



Crete catalogue - All depths - All magnitudes 




5 10 15 

log(x) 



FIG. 9: Weibull plot of CSS interevent times. The straight line indicates perfect Weibull 
dependence; m est := rh is the maximum likelihood estimate of the Weibull modulus from the 

data. 



B. Analysis of CSS interevent times for events above Ml, c 

Below, we focus on interevent times between events above thresholds that exceed the magni- 
tude of completeness. The resulting sequences are considered to be complete since they do not 
suffer from resolution issues. We show the Weibull plots for different magnitude thresholds in 



Fig. 10a The $(r; M > Ml jC ) curves show deviations from the Weibull in both the lower and 
upper tails. The curvature of the upper tail appears to reverse its sign as M L:C changes from 2.7 
to 2.9. Nevertheless, the Weibull dependence remains a useful first approximation. In particular, 
for all the thresholds considered, the Kolmogorov - Smirnov test does not reject the null hypoth- 
esis that the CDF is the Weibull with the respective best-fit (maximum likelihood) parameters at 
significance level 5%. 

In Table|n]we list the Weibull ITD statistics using different magnitude thresholds. Based on the 
confidence intervals for rh in Table |Il| there is no significant variation of rh with Ml )C . In contrast, 



r s increases with M L c as intuitively expected. In Fig. 10b we plot log(f s ) versus M L c G [2.3, 3.3] 
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(b) CSS time scales 



FIG. 10: (a) Weibull plot of CSS interevent times for different magnitude thresholds, (b) Semilog 



plot of empirical (circles) and theoretical -according to Eq. ([18])- (broken line) Weibull time scale 

versus Ml c - 



The best-fit (minimum least squares) line passing through the data is log(f s ) = A — pu Ml, c with 
A = 5.47, pm = 2.15, while the 95% confidence interval is [1.76, 2.54] . The value of pM predicted 



by the analysis in IV B is p* M = b log(10) « 1.96 in good agreement with the experimental p M . 



VI. CONCLUSIONS 



We propose a stochastic stick - slip model for the earthquake interevent times distribution that is 
based on the evolution of shear stress in the Earth's crust and the assumption that the crust strength 
follows the Weibull distribution, as is common for brittle materials. The current model differs from 
statistical approaches that test various empirical distributions, from universal scaling laws based on 
the concept of criticality ir7T Ji~0~ll23ll . from approaches rooted in probability theory BTi . and from 
stochastic branching models [fTTTl that incorporate by construction the Gutenberg - Richter and 
Omori laws. Our model is applicable to tectonic earthquakes that are responsible for the majority 
of the global earthquake activity. The interevent times distribution depends on both the crustal 
strength and the stress accumulation models. We show that the commonly used Weibull (including 
the exponential) and log-Weibull interevent times distributions are derived from the stochastic stick 
- slip model using, respectively, power-law (including linear) and logarithmic dependence of the 
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TABLE II: Number of events (N>) with magnitude above M L)C (2.3 - 3.3 M L ), ITD Weibull 
parameters rh and f s and their 95% confidence intervals, [mi, m 2 ] and [r s ,i, T s p] respectively. 
Weibull modulus is rounded to the second decimal place, while f s is rounded to the third most 

important digit. 



M L , C 


7V> 


rh 


t s (sec) 


[mi,m 2 ] 


[t~s,1, T a> 2] 


2.3 


629 


0.67 


3.46 10 4 


[0.66, 0.74] 


[3.07 10 4 , 3.89 10 4 ] 


2.5 


415 


0.70 


5.35 10 4 


[0.65, 0.75] 


[4.63 10 4 , 6.18 10 4 ] 


2.7 


335 


0.72 


6.80 10 4 


[0.66, 0.78] 


[5.81 10 4 ,7.9710 4 ] 


2.9 


177 


0.71 


1.25 10 5 


[0.63, 0.80] 


[1.04 10 5 , 1.62 10 5 ] 


3.1 


140 


0.70 


1.63 10 5 


[0.62, 0.80] 


[1.27 10 5 , 2.10 10 5 ] 


3.3 


70 


0.67 


3.14 10 5 


[0.55, 0.81] 


[2.17 10 5 , 4.56 10 5 ] 



stress accumulation function on time. We focus on power-law stress accumulation, which leads 
from Weibull strength statistics to Weibull interevent times. We also demonstrate that deviations 
from the Weibull dependence arise due to limited resolution (i.e., sampling of events exceeding a 
finite magnitude threshold), and random fluctuations of the accumulation rate. 



We also derive the scaling relation ( [18] ) that links the magnitude threshold with the scale of 
the respective Weibull interevent times distribution. This relation is based on the assumption that 
the Weibull model is an acceptable approximation of the ITD at finite magnitude thresholds (i.e., 
neglecting the lower-tail deviations) and on the Gutenberg - Richter law. The utility of the scaling 
relation is greater if the Weibull modulus can be assumed to remain relatively constant as the 
magnitude threshold changes. Then, this relation can be used to infer the Weibull time scales for 
larger-magnitude events based on the time scales of smaller- magnitude, more frequent, events. 

In future work we will investigate dynamic stick - slip models with controlled stress accumu- 
lation rates that will allow simulating seismic events. We will also focus on the deviations from 
the Weibull distribution in the upper tail. Another direction that will be pursued is the interplay 
between empirical scaling laws of fault and earthquake parameters (e.g., fault size distribution, 
seismic stress drop versus fault size and earthquake magnitude) and the ITD of a fault system. In 
particular, a model for the behavior of a fault system should involve an average over the sizes and 
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characteristic scales of the faults it comprises. Such an average should account for the empirical 
statistical facts that govern the fault population. While we believe that the interplay between the 
fault strength and the stress accumulation function is a universal factor controlling the ITD, the 
specific stress accumulation scenarios proposed and investigated herein do not exhaust the pos- 
sible functional forms of stress accumulation. In particular, seismic events driven by invasion of 
fluids in a fault system may be generated by more local and erratic stress accumulation patterns. 
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Appendix A: Effective strength distribution for linear depth dependence of Weibull scale parameter 



We assume that the strength scale depends on the depth as S s (h) = S + qh. Then, us- 
ing the variable transformation z = (So + qh)~ ms , and dz = —m s q(So + qh)~( m3+1 ^ dh = 
—m s q Z 1+1 / ms dh, i.e., dh = —^- q z~( 1+1 / ms \ the integral representing the survival function in 
Eq. |7]) is expressed as follows 



R* s (s) = -L_ — fdi^e-^, 

sv m a g(/i 2 -/n) J Z2 



where Zi = (So + qhj) ms , i — 1, 2. We now introduce the following definitions: 25h = h 2 — hi, 
S = Sl ~^ 32 , and the dimensionless variables: s = S/S, X s = q5h/S, u = zS m % and u lj2 = 

/ \ "iri s 

( j-r^- J , in terms of which the integral takes the form of Eq. ([8]). 
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1. Evaluation of -Rg(s) for m s > 1 



We evaluate the integral in Eq. ([8]) using integration by parts, using the survival function 



R*(s) = 1 - F*(s). This leads to 



R* s (s) 



-1 
2Al 



2 A, 



du (u- 1/ms )' e~ bu 
(u- 1/ma e~ bu )\ Ul +b I 1 duu- 1/ms e~ bu 



The term that originates from the boundary points (recalling that b = s m3 ) is equal to 



;i + a 5 

2 A, 



• 3 \ m S 
1 + — 



1-A 



2A< 



si - — V- 



For the remaining integral, we use the variable transformation z = bu, in terms of which we obtain 



« i 



h I duu- 1,ms e- bu = b 1/ms / dzz A '-- ■ 

' U2 J Z2 

b v™. [7(1 _ 1/ms> Zi) _ 7(1 _ 1/ms) Z2)] 



where z\ 2 



and j(a, x) is the lower incomplete gamma function defined by 7(0, x) 



duu a 1 e Hence, collecting the relevant terms we obtain 



R* s (s) 



;i + a a 

2A S 



. 1+ S A S ) 



2 A, 



s 

2X 



7 1 



-7 1 



m s ' (1- \ s ) m * 



m s (1 + A, 



(Al) 



The absolute difference between the CDF that is based on the numerical integration of Eq. ([8]) and 
that obtained from the explicit solution ( |A1[ ) is less than 1.4 x 1CT 8 over the parameter space of 
Fig.[TJ 



2. Evaluation of Ri,(s) for m s < 1 

For m s < 1 partial integration is not as efficient as for m a > 1, because the incomplete gamma 
function is not defined for 1 — ^- < 0; an extension of the gamma function is possible by means of 
partial integration, which reduces the incomplete gamma function of a negative argument to one 
with a positive argument. Nevertheless, for m s « 1 it follows that 1 — << —1, and repeated 



30 




FIG. 11: Dependence of the terms a n in Eq. ( |A3[ ), on n = 0, 1, ... 15 and s = 0, 0.5, ... 10, for 

parameter values \ s = 0.2, and m s = 0.7. 



partial integrations are required. Alternatively, a series expansion of the integral in Eq. ([8]) is 
obtained by means of the Taylor expansion of the exponential function er hu = X]^Lo(~ l)"^^ - 
(where b = s ms ). The survival function is then given by 

I x cx \ ( l) n g nm 



1 

2~X 



2m, A. 



E 

n=0 



n! 



duu n - {1+1/ms) 



E 

n=0 



(-l) n s nm ' (1 - A s ) 



1— n m s 



-d+A 4 



v 1— n m 3 



111 



nm s 



(A2) 



The above alternating series can be expressed as -Rg(s) = Yl^Lo(~^) n a «' where 



(1-A.) 



1— n m a 



'i + K 



, 1— n m s 



nl 



(A3) 



2 A s m s — 1) 

The sequence a n has a maximum for an integer n* which depends on s, A s , m s , as shown graph- 
ically in Fig.[TIJ Fixing these three parameters, we can write R%,(s) = Rg\s) + R^ n \s) where 

rV\s) is the finite series R ( s f \s) = Enlo(- 1 )" a « and R t\ s ) = EZii^T a 'n, 

where a' n = a n * +n . For n > 1 it holds that a' n < a' n+1 . Based on the alternating series test, 
R^ n \s) converges. In addition, if Rg n \s) is truncated after M terms, the absolute value of the 
remainder is less than aju+i- The absolute difference between the CDF that is based on the nu- 
merical integration of Eq. ([8]) and that obtained from the series ( |A2[ ) truncated at n = 100 is less 
than 2.2 x 10~ 9 over the parameter space of Fig. [T| 
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